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In  this  paper,  MAPLE  18  codes  are  used  to  utilize  Variational 
Iteration  method  for  the  numerical  solution  of  two  species  Lotka- 
Volterra  prey-predator  interaction  species  which  are  governed  by  a 
system  of  nonlinear  differential  equations.  Two  examples  are 
provided  to  show  the  ability  and  reliability  of  the  method.  The 
obtained  approximate  solution  shows  that  Variational  Iteration 
Method  (VIM)  is  powerful  numerical  technique  for  solving  a  system 
of  nonlinear  differential  equation,  which  can  be  easily  applied  to 
other  nonlinear  problems  in  biomathematics.  This  technique  has 
shown  to  be  very  effective  and  yields  accurate  results. 
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1.  Introduction 

System  of  differential  equation  has  a  wide  field  in  pure  and  applied  mathematics  such  as 
transportation  problem,  economics  mathematics,  meteorology,  biomathematics  and  engineering.  All 
of  these  disciplines  are  concerned  with  the  properties  of  differential  equations  of  various  types  which 
emphasizes  the  rigorous  justification  and  interpretation  of  natural  phenomena.  It  plays  an  important 
role  in  modelling  virtually  every  physical,  technical,  or  biological  process,  from  celestial  motion,  to 
bridge  design,  to  interactions  between  neurons.  Nonlinear  differential  equations  are  used  for 
describing  many  phenomena  in  the  real  world  as  prey  predator  interactions.  Prey  predator  models 
are  classified  as  one  of  the  most  important  applications  in  applied  mathematics,  thus  many  numerical 
and  semi-analytical  methods  are  developed  for  finding  the  solution  of  these  problems  by  many 
researchers  [1,  2,  3]. 

Relevant  biological  models  may  involve  interactions  (from  the  biochemical  to  the  ecosystem  level). 
One  of  the  first  interactions  model  in  population  dynamics  was  introduced  in  the  beginning  of  the 
20th  century  by  Alfred  Lotka,  an  American  biophysicist  (1925)  and  Vito  Volterra,  an  Italian 
mathematician  (1926)  models  interactions  between  preys  and  predators.  The  population  sizes  are 
denoted  by  x(t)  preys  and  y(t)  predators  at  timet  >  0.  He  assumed  that  population  change  of  one 
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specie  depends  on  its  current  population,  reproduction  rate  and  interactions  with  other  species 
(predator  and  prey)  [4]. 

In  recent  years,  research  on  Lotka-Volttera  system  of  nonlinear  differential  equations  have  been 
extensively  investigated  in  the  literature.  [5]  proposed  a  new  method  for  the  explicit  integration  of 
Lotka-Volterra  equations,  the  authors  in  [6]  solved  stochastic  Lotka-Volterra  equations  via 
operational  matrices,  a  study  on  competitive  Lotka-Volterra  model  in  random  environments  was 
examined  and  studied  by  [7].  [8]  studied  and  observed  the  population  dynamics  of  a  three-species 
Lotka-Volterra  model  of  two  predators  and  their  prey,  [9]  applied  differential  transformation  method 
for  solving  prey  predator  model  with  boiling  type  I.  [10]  studied  a  series  solutions  of  two  species 
Lotka-Volterra  Equations  by  Adomian  Decomposition  and  Homotopy  Perturbation  Methods  and 
[11]  studied  Predator-Prey  system  in  a  polluted  environment. 


1.1  Lotka-Voltarra  Model 

Lotka-Volterra  model  also  known  as  the  predator-prey  equations,  in  deterministic  subclasses,  are 
well-known  and  have  been  an  active  area  of  research  concerning  ecological  population  modeling 
and  economic  modeling.  These  type  of  equations  is  so  attractive  in  the  terms  of  population  dynamics 
of  species  competing  (conflict)  and  the  logistic  population  model.  Thus,  the  Lotka-Volterra  model 
in  case  of  two  species  is  a  prey  predator  equation  which  is  defined  as  follows: 
f  dx 


—  =  ax{t)  -  /?x(t)y(t) 

<^y 

+  TX{t)y{t) 


(1) 


subject  to  initial  conditions 

xito)  =  A 
y(to)  =  B 

where  the  function  x(t)  represents  the  populations  of  prey  at  time  t,  and  also  the  function  y(t) 
represents  the  populations  of  predator  at  time  t.  All  of  the  parameters  a,  p,  o),  r  are  non-negative 
constants.  The  parameter  a  represents  the  per  capita  reduction  in  prey  per  predator.  The  parameter 
represents  death  rate  per  encounter  of  prey  due  to  predation.  Moreover  if  p  is  the  only  decreasing 
factor  for  the  prey  population,  then  prey  will  be  eaten  by  predators.  The  parameter  o)  represents  the 
per  capita  increase  in  predator  per  prey,  moreover  if  o)  is  the  only  increasing  factor  for  the  predator 
population,  then  the  population  growth  is  proportional  to  the  food  available.  The  parameter  t 
represents  mortality  rate  of  predator  and  A,  B  are  constants. 

In  reality,  if  the  prey  population  is  large,  the  predators  will  have  more  food  to  support  a  larger 
population.  However,  when  the  predator  population  grows  too  large,  the  prey  begins  to  die  off.  This 
will  result  in  a  decrease  in  the  predators.  This  trend  continues  as  time  goes  on,  implying  a 
stable  coexistence  of  the  two  populations. 

The  main  objectives  of  this  paper  is  to  present  and  employ  a  MAPLE  18  software  codes  for 
variational  iteration  method  proposed  in  [12]  and  to  overcome  the  mathematical  stress  of  integral 
involve  in  implementation  of  VIM. 

This  paper  is  organized  as  follows.  In  section  1,  we  introduced  briefly  the  Lotka-Volterra  model  and 
their  parameters  notation  related  to  Predator-prey  interaction  behaviors.  Moreover,  stability  and 
equilibrium  of  the  model  are  discussed.  In  section  2,  we  discussed  and  formulate  MAPLE  18 
Variation  iteration  scheme.  Section  3,  the  variational  iterative  method  was  applied  to  solve  predator- 
prey  Equation  (1)  while  in  section  4,  the  numerical  results  table  and  graphs  are  reported.  Einally, 
section  5  provides  the  discussion  and  conclusion. 
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1.2  The  Dynamic  Behavior  of  the  Lotka-Volterra  Model 

One  of  the  main  properties  of  dynamic  systems  is  stability.  The  stability  is  studied  to  determine 
some  properties  of  solutions  or  system  of  differential  equations.  Consequently  the  dynamic  behavior 
of  the  model  will  be  discussed.  For  the  model  equilibrium  points,  we  set  the  right  hand  side  of  (1) 
to  zero. 


(a  -  py{t))x{t)  =  0 
(— (x)  +  Tx(t))y(t)  =  0 

to  obtain  x  =  0,y  =  0orx  =  — ,y  =  ^.  Then  system  (1)  has 

T  p 


i. 

ii. 


the  trivial  equilibrium,  £'o(xo,yo)  =  (0,0), 
the  non-tri vial  equilibrium  point,  F^Cx^.y*)  = 


(3) 


To  investigate  the  stability  of  each  equilibrium  point,  we  evaluate  the  Jacobian  matrix  of  (1)  at  each 
of  these  equilibrium.  The  Jacobian  matrix  of  system  (1)  is  given  by 

—  Py  —px 

ry  —0)  +  TXj 


J(x,y)  = 


(4) 


1.3  At  the  Trivial  Equilibrium 

We  obtain  from  (4)  that 

ASo)  =  (o  -„)• 

Thus  the  two  eigenvalues  of /(Eq)  are  =  a  and  A2  =  —o).  Therefore,  Eq  is  a  saddle  point  since 
one  of  the  eigenvalues  is  positive  and  the  other  is  negative. 

The  corresponding  eigenvectors  are  the  two  axes  x  and  y.  Then  using  the  standard  computation  of 
eigenvectors,  we  have  for  the  eigenvector  associated  with  =  a 

(0 

so  that  (a  +  a))y  =  0  and  so,  y  =  0  which  is  the  x  —axis.  Therefore,  the  x  component  of  a 
perturbation  away  from  but  still  near  the  Eq  will  grow  exponentially  at  a  rate  a. 

Similarly,  x  =  0,  which  is  the  y  —  axis,  is  the  eigenvector  associated  with  A2  =  —o).  Therefore, 
the  y  component  of  a  perturbation  away  from  but  still  near  the  Eq  will  shrink  exponentially  at  a  rate 

-CO. 


1.4  At  the  uou-trivial  equilibrium  (Coexisteuce  Equilibrium) 

The  Jacobian  matrix  in  (4)  becomes: 


(7) 


The  eigenvalues  of  J(E*)  areAi  2  =  +iVaw.  This  show  that  the  trace  of  this  matrix  is  zero  and  the 
determinant  isam  >  0.  Thus  E*  is  a  centre  and  is  balanced  at  the  knife-edge  between  stable  and 
unstable  oscillations.  Therefore  system  (1)  is  structurally  unstable  since  any  slight  change  to  the 
structure  of  the  equations,  especially  changing  the  nonlinear  terms,  could  tip  the  balance  between 
stability  or  instability,  depending  on  how  changes  to  the  structure  of  the  equations  affects  the  real 
part  of  the  eigenvalues.  In  a  structurally  unstable  system,  slight  modifications  to  the  form  of  the 
equations  alter  the  stability. 

The  counterclockwise  circulation  of  vectors  near  the  coexistence  equilibrium  E*  is  as  a  result  of  the 
purely  imaginary  eigenvalues  of  (E*).  Any  perturbation  near  this  equilibrium  causes  the  system  to 
oscillate  around  it  infinitely  in  a  closed  orbit,  neither  growing  away  from  the  equilibrium  nor 
returning  to  it.  The  larger  the  perturbation,  the  greater  the  amplitude  of  the  circulation,  see  Figures 
4  and  6  for  example. 
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2.  Description  of  the  Variation  Iteration  Method  (VIM)  and  Solution  Approach 

The  basic  idea  of  the  variational  iteration  method  is  to  construct  an  iteration  procedures  based  on 
correction  functional  that  include  a  generalized  Lagrange  multiplier  [13-14].  The  VIM  was 
proposed  where  the  value  of  the  multiplier  was  chosen  using  variational  theory  so  that  each 
improves  the  accuracy  of  the  solution  [15]. The  initial  approximation  i.e.  trial  function  usually 
includes  unknown  coefficient  which  can  be  determined  to  satisfy  any  boundary  and  initial 
conditions.  VIM  does  not  require  specific  transformation  for  nonlinear  terms  as  required  by 
other  techniques  and  is  now  widely  used  by  many  researchers  to  study  autonomous  ordinary 
differential  equation,  Integro-differential  systems,  Linear  Helmholtz  partial  differential 
equation  and  other  fields  [16-21].  In  this  method  the  solution  is  given  in  an  infinite  series  usually 
convergent  to  an  accurate  solution.  According  to  the  variational  iteration  method  we  consider  the 
following  general  differential  equation  of  the  form: 

Lp  +  Np  =  g(t)  (8) 

where  L  is  a  linear  operator  V  is  a  nonlinear  operator  and  g(t)  is  an  inhomogeneous  term.  We  can 
construct  a  correctional  function  as  follows 

Xn+i  =  Xn(t)  f^A{LXn(s)  +  N^(s)  -  g(s)}  ds  (9) 

Where  A  is  a  Lagrangian  multiplier  which  can  be  identified  optimally  via  variational  theory  [22]. 
The  subscript  n  denotes  the  nth  approximation  Consider  the  stationary  condition  of  the  above 
correction  functional  then  the  Lagrange  multiplier  can  be  expressed  as 

(10) 

Where  q  is  the  highest  order  of  the  differential  equation. 

2.1  Maple  18  Coded  Variational  Iteration  Method 

In  order  to  formulate  the  general  variational  iteration  approach  on  MAPLE  1 8  software,  we 
consider  Equations  (1)  and  (2)  and  develop  the  VIM  schemes  as  follows: 

MAPLE  18  Coded  Variational  Iteration  Scheme 
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r  restart: 

a  :=  R^)  P  ■■=  Rz',  0)  ■■=  Rz',  t  :=  R^-,  A  ■■=  (-1); 

^0  ==  A;  yo  •=  5;  Uq  :=  diff(xo,t);  %  :=  diff(yQ.t); 

:=  Xq  +  (A)  *  mt((ao  —  (a  *  Xq)  +  (^  *  Xq  *  yo))>  t  =  0  ...  t); 
Xii  :=  value(xi);  bo  ■=  di//(xii,  t); 
yi  ==  yo  +  W  *  mt{(ai  +  (o)  *  yo)  -  (t  *  Xo  *  yo)),  t  =  0  ...  t}; 
:=  vaZue(yi);  :=  diffiy-^^^.t); 

Xz  •■=  Xii  +  (A)  *  int{(6o  -  (a  *  x^)  +  (/?  *  x^  *  y^)),  t  =  0  ...  t}; 

Xi2  ==  value{xz)]  Cq  ■=  diffix-^z.t)-, 
y2  ==  yii  +  W  *  +  (co  *  yii)  -  (t  *  Xii  *  yii)),  t  =  0  ...  t}; 

yi2  ==  valueiyz);  ■=  di//(yi2,t); 

Xs  :=  Xi2  +  (A)  *  int{{co  -  (a  *  X12)  +  (/?  *  X12  *  yi2)),  t  =  0  ...  t}; 

Xi3  :=  value(xz);  do  •=  diff(xiz,t)', 
ys  ==  yi2  +  W  *  mt{(ci  +  (co  *  yi2)  -  (t  *  x^  *  yi2)),  t  =  0  ...  t}; 
yi3  :=  valueiyz)]  di  ==  diffiy^^.t); 


Xm  ■=  Xim-l  +  U)  *  int{{so  “  («  *  ^Ti^-i)  +  (/?  *  *  yim-l))- 1  =  0  ...  t}; 

^im  ==  ValueiXm)-,  Zo  ■=  diffix^rn>t)> 

Vm  •=  yim-1  +  W  *  mt{(si  +  (o)  *  yi^_i)  -  (t  *  Xi^_i  *  yim-i))- 1  =  0  ...  t]; 
yim  ■=  valueiy^y,  ■■=  diffiy^^.t)-, 

X(t)  :=  evalf  {x-im)) 

I  y(t)  :=  evalfiy^rn)', 

where  m  is  the  number  of  iterations,  R^,  Rz,  R3,  R4  are  parameter  constants  and  A,  B  are  constants. 
Thus,  x(t)and  y(t)  are  the  series  solutions  at  mth  iteration. 


3.  Numerical  Experiment 


In  this  section,  we  examine  the  predator  and  prey  growth/decay  of  the  two  species  (self-interaction) 
as  well  as  their  interaction.  We  investigate  bahaviours  of  the  parameters  a,  p,  co  and  x  which  play  a 
key  role  in  determining  the  system  behaviours.  Two  examples  are  considered  subject  to  specific 
initial  condition  of  the  model.  Using  the  MAPLE  18  coded  Variation  Iteration  scheme,  we  obtained 
the  numerical  solution  for  Equation  (1)  and  initial  condition  (2)  at  iteration  step  m  =  10  as  shown 
in  Tables  1  to  6. 


Example  1  Suppose  a  is  the  per  capita  reduction  in  prey  per  predator  and  P  is  the  death  rate  per 
encounter  of  prey  due  to  predation  [23]. 


a  =  3.029850 
/?  =  3.029850 
0)  =  1.967217 
T  =  2.295942 
A  =  -1.00000 
Xo  =  1.187100 


Substitute  the  above  parameters  into  algorithm  (1 1),  we  have  the  following  solutions: 


when  {a  <  P) 

^  when  (a  >  P)  ^ 

a  =  3.029850 

a  =  4.094132 

p  =  4.094132 
m  =  1.967217 

T  =  2.295942 

P  =  3.029850 
m  =  1.967217 

T  =  2.295942 

A  =  -1.00000 

A  =  -1.00000 

Xo  =  1.187100 

Xo  =  1.187100 

.yo  =  0.740047 

.yo  =  0.740047  . 
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Table  1;  (a  <  P) 


t 

Analytical 
Solution  x(t) 

VIM 

Solution  x(t) 

Analytical 
Solution  y(t) 

VIM 

Solution  y(t) 

0.0 

1.187100000000 

1.187100000000 

0.7400470000000 

0.7400470000000 

0.1 

1.173293314174 

1.173293331000 

0.7975048380086 

0.7975047610000 

0.2 

1.132681174388 

1.132681233000 

0.8540331977419 

0.8540330912000 

0.3 

1.069749554336 

1.069749539000 

0.9036493434085 

0.9036492527000 

0.4 

0.992376686741 

0.992376610300 

0.9407098212139 

0.9407097698000 

0.5 

0.909715706923 

0.909715551500 

0.9612974253317 

0.9612976897000 

0.6 

0.829963573965 

0.829961925000 

0.9640576734517 

0.9640598752000 

0.7 

0.758989599935 

0.758979098800 

0.9501583685533 

0.9501676352000 

0.8 

0.700090581059 

0.700045711900 

0.9225636508131 

0.9225907487000 

0.9 

0.654510497041 

0.654363278700 

0.8850829288848 

0.8851423567000 

1.0 

0.622205729269 

0.621805168400 

0.8415741110666 

0.8416743875000 

Table  2:  (a  >  jg) 


t 

Analytical 
Solution  x(t) 

VIM 

Solution  x(t) 

Analytical 
Solution  y(t) 

VIM 

Solution  y(t) 

0.0 

1.187100000000 

1.187100000000 

0.7400470000000 

0.740047000000 

0.1 

1.413308431481 

1.413308491000 

0.8192192423347 

0.819219212200 

0.2 

1.629422801606 

1.629422522000 

0.9548299457312 

0.954830107800 

0.3 

1.783965142331 

1.783964837000 

1.1626432476204 

1.162643502000 

0.4 

1.812364112012 

1.812363259000 

1.4474459753223 

1.447445942000 

0.5 

1.674587231402 

1.674575344000 

2.0850176428230 

2.085017642720 

0.6 

1.401596690416 

1.401562584000 

2.0850176428292 

2.085017641320 

0.7 

1.085902388192 

1.086180815000 

2.2783956366840 

2.278395618000 

0.8 

0.811351279883 

0.811351460690 

2.3239540305871 

2.323951416000 

0.9 

0.610020458917 

0.610020361200 

2.2440744293156 

2.244074124000 

1.0 

0.476079350030 

0.476079236100 

2.0857367907384 

2.085736623000 

Table  3:  (a  =  P) 


t 

Analytical 
Solution  x(t) 

VIM 

Solution  x(t) 

Analytical 
Solution  y(t) 

VIM 

Solution  y(t) 

0.0 

1.187100000000 

1.187100000000 

0.740047000000 

0.740047000000 

0.1 

1.272187665251 

1.272187275000 

0.806474171253 

0.806474407600 

0.2 

1.331974163143 

1.331973873000 

0.893823428997 

0.893823624800 

0.3 

1.354102768524 

1.354102583000 

1.000208677394 

1.000208771000 

0.4 

1.330095316143 

1.330095126000 

1.119123022586 

1.119122920000 

0.5 

1.259716623148 

1.259716980000 

1.238564988790 

1.238564225000 

0.6 

1.152840256063 

1.152845461000 

1.342772096989 

1.342768205000 

0.7 

1.026594271894 

1.026627603000 

1.416793455406 

1.416793581000 

0.8 

0.899120517439 

0.899244376100 

1.451520906161 

1.451521131000 

0.9 

0.784059292436 

0.784362541500 

1.445978046355 

1.445913412000 

1.0 

0.688601146327 

0.688090398100 

1.405950208844 

1.405941241000 
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Example  2  Suppose  a  is  the  per  capita  reduction  in  prey  per  predator  and  ^  is  the  death  rate  per 
encounter  of  prey  due  to  predation  [24]. 


when  (a  > 

when  (a  < 

when  (a  =  P) 

a  =  0.100 

a  =  0.0014 

a  =  0.100 

/?  =  0.0014 

p  =  0.100 

p  =  0.100 

o)  =  0.080 

u)  =  0.080 

0)  =  0.080 

T  =  0.0012 

T  =  0.0012 

T  =  0.0012 

A  =  -1.000 

A  =  -1.000 

A  =  -1.000 

Xq  =  4.000 

Xq  =  4.000 

Xq  =  4.000 

.  yo  =  9.000 

.  yo  =  9.000 

.  yo  =  9.000 

Substitute  the  above  parameters  into  algorithm  (1 1),  we  have  the  following  solutions: 
Table  4:  (a  >  P) _ _ _ _ 


t 

Analytical 
Solution  x(t) 

VIM 

Solution  x(t) 

Analytical 
Solution  y(t) 

VIM 

Solution  y(t) 

0.0 

4.000000000000 

4.000000000000 

9.000000000000 

9.000000000000 

0.1 

4.035132287030 

4.035132286000 

8.932592638256 

8.932592638000 

0.2 

4.070611403006 

4.070611402000 

8.865727698739 

8.865727700000 

0.3 

4.106440756278 

4.106440756000 

8.799400925152 

8.799400926000 

0.4 

4.142623789899 

4.142623790000 

8.733608094714 

8.733608095000 

0.5 

4.179163981615 

4.179163982000 

8.668345018151 

8.668345018000 

0.6 

4.216064844061 

4.216064844000 

8.603607539534 

8.603607540000 

0.7 

4.253329953021 

4.253329924000 

8.539391509681 

8.539391537000 

0.8 

4.290962896869 

4.290962804000 

8.475692832270 

8.475692921000 

0.9 

4.328967276114 

4.328967104000 

8.412507466043 

8.412507636000 

1.0 

4.367346731478 

4.367346480000 

8.349831399216 

8.349831655000 

Table  5:  (a  <  P) 


t 

Analytical 
Solution  x(t) 

VIM 

Solution  x(t) 

Analytical 
Solution  y(t) 

VIM 

Solution  y(t) 

0.0 

4.000000000000 

4.000000000000 

9.000000000000 

9.000000000000 

0.1 

3.657473020744 

3.657473258000 

8.932387296668 

8.932387294000 

0.2 

3.346536796070 

3.346536998000 

8.864935239948 

8.864935237000 

0.3 

3.064097573818 

3.064097837000 

8.797679574554 

8.797679571000 

0.4 

2.807379952412 

2.807380177000 

8.730651950326 

8.730651947000 

0.5 

2.573892073952 

2.573892359000 

8.663880372174 

8.663880368000 

0.6 

2.361396256269 

2.361396506000 

8.597389579842 

8.597389576000 

0.7 

2.167881361179 

2.167881643000 

8.531201405234 

8.531201402000 

0.8 

1.991539481074 

1.991539751000 

8.465335073997 

8.465335070000 

0.9 

1.830744135389 

1.830744404000 

8.399807487591 

8.399807484000 

1.0 

1.684031442829 

1.684031744000 

8.334633466917 

8.334633463000 

Table  6:  (a  =  jg) 


t 

Analytical 
Solution  x(t) 

VIM 

Solution  x(t) 

Analytical 
Solution  y(t) 

VIM 

Solution  y(t) 

0.0 

4.000000000000 

4.000000000000 

9.000000000000 

9.000000000000 

0.1 

3.693713901517 

3.693714070000 

8.932407277368 

8.932407275000 
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0.2 

3.413183652004 

3.413183941000 

8.865010279466 

8.865010275000 

0.3 

3.156082329242 

3.156082424000 

8.797838192760 

8.797838190000 

0.4 

2.920304355828 

2.920304706000 

8.730917063216 

8.730917058000 

0.5 

2.703946085963 

2.703946207000 

8.664270073277 

8.664270071000 

0.6 

2.505282320509 

2.505282688000 

8.597917875002 

8.597917871000 

0.7 

2.322752234866 

2.322752401000 

8.531878791541 

8.531878789000 

0.8 

2.154939736038 

2.154940050000 

8.466169095953 

8.466169091000 

0.9 

2.000562151892 

2.000562388000 

8.400803173218 

8.400803170000 

1.0 

1.858455074810 

1.858455294000 

8.335793735914 

8.335793731000 

4.  Graphs  Representation 


I - xm - y(t)1 

Figure  1  Numerical  solutions  of  the  Predator  y(t)  and  Prey  x(t)  when  a  <  P  Example  1 


Figure  2  Unsteady  states  of  phase  planes  for  Predator-Prey  when  a  <  P  Example  1 


Prey  y('*)  Fneda^Qrl 


Figure  3  Numerical  solutions  of  the  Predator  y(t)  and  Prey  x(t)  when  a  >  P  Example  1 
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Figure  4  Steady  states  of  phase  planes  for  Predator-Prey  when  a  >  P  Example  1 


I - x(t)  Prey - y(t)  Predator| 


Figure  5  Numerical  solutions  of  the  Predator  y(t)  and  Prey  x(t)  when  a  =  P  Example  1 


Figure  7  Numerical  solutions  of  the  Predator  y(t)  and  Prey  x(t)  when  a>  P  Example  2 
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Figure  8  Unmutualistic  interactions  phase  for  Predator-Prey  when  a>  P  Example  2 


I  Prey  —  y<*j  Fre<3ator| 


Figure  9  Numerical  solutions  of  the  Predator  y(t)  and  Prey  x(t)  when  a  <  P  Example  2 


Figure  10  Mutualistic  interactions  phase  for  Predator-Prey  when  a  <  P  Example  2 
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Figure  11  Numerical  solutions  of  the  Predator  y(t)  and  Prey  x(t)  when  a  =  ^  Example  2 


Figure  12  Unmutualistic  interactions  phase  for  Predator-Prey  when  a  =  P  Example  2 


3.1  Discussion 

The  numerical  solutions  for  six  cases  of  interaction  behaviors  between  Predator-Prey  are  presented 
in  Tables  1  to  6  while  Figures  1  to  6  show  plots  relationship  of  Predator-Prey  for  example  one  and 
Figures  7  to  12  depict  the  numerical  solution  and  interaction  behaviors  of  Predator-Prey  for  example 
two.  Finally,  the  study  have  revealed  the  interactions  behaviors  between  two  species  (steady, 
unsteady,  mutualistic  and  unmutualistic  interactions)  from  computational  analysis. 

4.  Conclusion 

This  article  highlights  the  feasibility  and  capability  of  MAPFE  18  software  codes  to  solve  nonlinear 
system  of  Lotka-Volterra  two  species  of  Predator-Prey  model.  The  conventional  variation  iteration 
method  was  employed  using  MAPLE  18  software  commands  to  overcome  the  rigorous 
computational  work  and  simplification  of  integrals  involves  during  iteration  process.  Two  examples 
are  consider  to  test  the  efficiency  of  the  algorithm  and  the  numerical  results  obtained  were  compared 
with  analytical  solutions  with  little  relative  errors.  The  advantage  of  the  MCVIM  over  the 
conventional  approach  is  faster,  easy  and  it  provides  an  efficient  numerical  solutions.  Finally,  the 
proposed  scheme  is  easy  to  implement  and  shows  a  good  agreement  with  analytic  results. 
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